String Method for the Study of Rare Events 



O 

o 



in 



> 

in 
in 
o 

(N 
O 



i 

c 

o 
o 



13 



Weinan E 1 , Weiqing Ren 2 , and Eric Vanden-Eijnden 2 
1 Department of Mathematics and PACM, Princeton University, Princeton, NJ 08544 
2 C our ant Institute, New York University, New York, NY 10012 

We present a new and efficient method for computing the transition pathways, free energy barriers, 
and transition rates in complex systems with relatively smooth energy landscapes. The method 
proceeds by evolving strings, i.e. smooth curves with intrinsic parametrization whose dynamics takes 
them to the most probable transition path between two metastable regions in the configuration space. 
Free energy barriers and transition rates can then be determined by standard umbrella sampling 
technique around the string. Applications to Lennard- Jones cluster rearrangement and thermally 
induced switching of a magnetic film are presented. 
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The dynamics of complex systems are often driven by 
rare but important events (for a review see e.g. Q ) . Well- 
known examples include nucleation events during phase 
transition, conformational changes in macromolcculcs, 
and chemical reactions. The long time scale associated 
with these rare events is a consequence of the dispar- 
ity between the effective thermal energy and typical en- 
ergy barrier of the systems. The dynamics proceeds by 
long waiting periods around metastable states followed 
by sudden jumps from one state to another. 

Sophisticated numerical techniques have been devel- 
oped to find the transition pathways and transition rates 
between metastable states in complex systems for which 
the mechanism of transition is not known beforehand 
P, 0, pj . With the exception of the transition path sam- 
pling technique j3] , most of these methods seem to require 
that the energy landscape be relatively smooth. One typ- 
ical example of such techniques is the nudged elastic band 
(NEB) method [Q. NEB connects the initial and the fi- 
nal states by a chain of states. The states move in a force 
field which is the combination of the normal component 
of the potential force and the tangential component of 
the spring force connecting the states. The spring force 
helps to evenly space the states along the chain. 

In this paper we propose an alternative approach for 
computing transition pathways, free energy barriers, and 
transition rates. We sample the configuration space with 
strings, i.e. smooth curves with intrinsic parametrization 
such as arclength, or energy weighted arclength which 
connect two metastable states (or regions), A and B. 
The string satisfies a differential equation which by con- 
struction guarantees that it evolves to the most probable 
transition pathway connecting A and B. One can then 
perform an umbrella sampling of the equilibrium distri- 
bution of the system in the hyperplanes normal to the 
string and thereby determine free energy barriers and 
transition rates. 

Consider the example of a system modeled by 

79=-W(g)+£(t). (1) 
where 7 is the friction coefficient, £(i) is a white-noise 



with (£j -(f)£fc(0)) = 27fesT5jfc5(t). The metastable states 
are localized around the minima of the potential V(q). 
Assuming V(q) has at least two minima, A and B, we 
look for the minimal energy paths (MEPs) connecting 
these states. By definition, a MEP is a smooth curve ip* 
connecting A and B which satisfies 



(VV) x (<p*) = 0, 



(2) 



where (VV) is the component of W normal to <p* . The 
MEPs are the most probable transition pathways for (|l|) 
since with exponentially high probability it is by these 
paths that the system switches back and forth between 
the states A and B under the action of a small thermal 
noise It is interesting to note that the solutions of 
(^|) also provide relevant information about the Langevin 
equation 



q = p, 

V=-W{q)- 1P + at)- 



(3) 



Indeed, the metastable regions for ([!]) and (||) coincide, 
and the transition pathways for ([|) can be easily deter- 
mined from the transition pathways for (^) because they 
traverse the same sequence of critical points. As a result 
the transition rates for (j|) for an arbitrary friction coef- 
ficient 7 can be obtained by considering the high friction 
evolution equation ([!]) - see (jl^) below. 

Let p> be a string (but not necessarily a MEP) connect- 
ing A and B. A simple method to find the MEP is to 
evolve ip according to 



= -(W) ± (<p), 



(4) 



where u 1 - denotes the normal velocity of <p, since station- 
ary solutions of (Q) satisfy (|^). For numerical purposes 
it is convenient to have a parametrized version of (Q), 
keeping in mind however that the parametrization can 
be arbitrarily chosen since both (||) and (||) are intrin- 
sic. Denote by <p(a,t) the instantaneous position of the 
string, where a is some suitable parametrization. Then 
we can rewrite (El) as 



< ft = -{VV{p>)) 1 - +rt, 



(5) 
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FIG. 1: Top figure: A transition pathway by which the cen- 
tral atom migrate to the surface in a seven-atom hexagonal 
Lennard Jones cluster in the plane. The pictures show succes- 
sive configuration corresponding to local minima of potential 
energy along the path. Bottom figure: the potential energy 
along the path in natural units. The full line corresponds to 
a simulation with N = 200 discretization points along the 
string, and the dashed line, N = 20. 

where for convenience we rcnormalizcd time i/7 t, 
(VV) = W — (VV-r)t, and t is the unit tangent vector 
along ip, t = fa/\(p a \. The scalar field r = r(a,t) is a 
Lagrange multiplier uniquely determined by the choice of 
parametrization. The simplest example is to parametrize 
ip by arclength normalized so that a = at A, a = 1 at 
B. Then (^) must be supplemented by the constraint 

(K|)a = 0. (6) 

which determines r Q. Other parametrizations can be 
straightforwardly implemented by modifying the con- 
straint (^|). For instance, a parametrization by en- 
ergy weighted arclength which increases resolution at 
the transition states is achieved using the constraint 
(f(V(<P))\<Pa\)a = 0, where f(z) is some suitable mon- 
itor function satisfying f'(z) < 0. In addition, the end 
points of the string need not be fixed and other boundary 
conditions can be used. 

Because of the intrinsic description of the string, it is 
very simple to implement an efficient algorithm which 
solves (||) using a time-splitting type of scheme. The 
string is discretized into a number of points which move 
according to the first term, — (VV^^))- 1 , at the right 
hand-side of ([5]). After a number of steps depending on 
the accuracy for the constraint (|^), a reparametrization 
step is applied to reinforce (||). This costs O(N) oper- 
ations where N is the number of discretization points 
along the string. At the reparametrization step it is also 




200 400 600 800 1000 1200 1400 1600 1800 2000 
number of iterations 



FIG. 2: The convergence history of the steepest descent 
method [(a)] and a limited memory version of Broyden's 
method [(b)] applied to the seven-atom cluster problem with 
N = 200 discretization points. We use maxo< a <i | ( W) | to 
measure accuracy. 

convenient to change N according to the accuracy re- 
quirement for the representation of the string. 

The method certainly bears some similarities with 
NEB since one can think of the introduction of the 
spring force in NEB as a way of ensuring equal-distance 
parametrization by a penalty method - NEB gives an 
evolution equation which, in the continuum limit, is sim- 
ilar to (^) but with r given by r = K(p aa • i where k is 
the artificial spring constant. As in other penalty meth- 
ods, this numerical procedure introduces stiffness into 
the problem if the penalization parameter, here the elas- 
tic constant «, is large and this limits the size of the 
time step. By using an intrinsic description, we elimi- 
nate this problem and speed up convergence. Further- 
more we gain the ability of using other parametrizations 
in a simple and flexible way. Finally, there is no simple 
way to change the number of discretization points along 
the chain in NEB. 

It is natural to ask how the string method compares to 
NEB in terms of performance. However such a compar- 
ison does not seem straightforward since it depends on 
the criteria. This is discussed in detail in ||. 

As a first example we look at the dynamics of seven 
atoms interacting via Lennard- Jones potential on the 
plane. This example has been studied in detail in [Q. 
In equilibrium the seven atoms form an hexagon. We are 
interested in the process in which the atom at the cen- 
ter migrates to an external position. The MEPs are not 
unique for this problem. One example of MEP obtained 
via the string method is shown in figure Q the critical 
points along this path coincide with the ones obtained in 
JtJ by transition path sampling. 

The string equation essentially amounts to finding 
the MEPs by the method of steepest descent except that 
we are not working with any explicit object function to 
minimize. We can use advanced numerical techniques for 
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solving nonlinear equations |J to accelerate convergence 
to the MEP. We have developed a limited memory version 
of Broyden's method where (||) is replaced by 



<pt = -G ± (Vy(v?))- L +rt. 



(7) 



Here G 1 - is a matrix determined on the fly to approxi- 
mate the inverse of the Hessian in the perpendicular hy- 
perplane; the approximation is based on the past history 
of ip and does not require to actually compute the Hes- 
sian (for details, see Q). In figure |2|, we compare the 
convergence history of the steepest descent method and 
the Broyden-accelerated method applied to the seven- 
atom cluster problem. The Broyden-accelerated method 
converges much faster. 

Once the MEP <p*(a) has been determined using 
Broyden-accelerated string method, free energy barriers 
and transition rates can be computed by standard um- 
brella sampling of the equilibrium distribution of the sys- 
tem in S* (a) , the hyperplane normal to ip* (a) . Consider 
first the free energy difference along the string defined as 
F(a) - F(0) = -k B Tln(Z(a)/Z(0)) where 



Z(a) = I e-P v ^d d q 

IS* {a) 



(8) 



is the partition function and (3 = 1/fceT. Using the 
identity J 9 In Z/dada — lnZ(a), we obtain from (|8|) 

F(a) - F(0) = f ((f* ■ W) ((t* ■ <p*) a , - i* a ■ q))da'. 
Jo 

(9) 

Here (■) is the ensemble average with respect to equilib- 
rium distribution restricted in the hyperplane S (a) and 
t* is the unit tangent vector along ip* . (fflis similar to 
the standard thermodynamic integration |10| but is bet- 
ter suited for numerical purposes. In practice, we use 
ergodicity and replace the ensemble average in (^|) by 
a time-average over the solution of an equation similar 
to (jl]) but restricted in the hyperplane S(a). Following 
Kramers original argument (see e.g. chap. 9 in [pj|), the 
transition rate can be expressed in terms of the free en- 
ergy as 



k A - 



2-\/A m | X s 



-f3AF 



tt(7+ vV +4|A fl |) 
where AF is free energy barrier along ip* 
AF = max (F(a) - F(0)) . 

0<Q<1 



(10) 



(11) 



and X m and X s are (inverse square of) characteristic time 
scales at the minimum and maximum of the free en- 
ergy along the transition path; A m and X s are given by 
\<p a \~ 2 F aa evaluated at a = and a = a s , respectively, 
where a s is the value at which the maximum in (|ll|) is 
attained |12fl. 



The transition rates along the MEP obtained earlier 
for the seven-atom cluster were evaluated by the string 
method, and are summarized in table B 

The string method can easily be generalized to infi- 
nite dimensional dynamical systems by introducing an 
appropriate norm in phase-space. As an example, we 
consider the problem of thermally induced switching of a 
magnetic film. This problem is of great current interest 
in the magnetic recording industry [ [15|. (For an intro- 
duction to micromagnetism, see e.g. p3, p^[ ; thermally 
induced switching is studied in Jj"iJ|). Landau-Lifshitz 
theory of micromagnetism provides an energy for a fer- 
romagnetic sample il which after suitable nondimension- 
alization, reads 

E[m]=A j \Vm\ 2 d 3 x+ / (f>(m)d 3 x + / \Vu\ 2 d 3 x, 
Jn Jn Jr 3 

(12) 

where m is the magnetization distribution normalized 
so that \m\ = 1. The three terms represent respec- 
tively energies due to exchange, anisotropy, and stray 
field. The potential u, defined everywhere in space, solves 
div (— Vtt + m) = 0, where m is extended as outside f2. 

Various switching pathways for ( |l2| ) were obtained us- 
ing the string method: two examples are shown in figure p| 
and the energy along these paths is shown in figure W. 
These paths illustrate two generic mechanisms for switch- 
ing in magnetic films. Path (a), which is more favorable 
in thin samples, proceeds by domain wall motion, inte- 
rior rotation, and switching of the edge domains. Path 
(b), which is more favorable for thicker films, proceeds 
by vortex nucleation, invasion of the sample and vortex 
expulsion. 

In conclusion, transition pathways and transition rates 
for complex systems with a relatively smooth energy 
landscape can be determined efficiently by evolving 
strings instead of points in configuration space. The in- 
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exact 


4.969 x KT 13 
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TABLE I: The rates for the various subprocesses in the tran- 
sition shown in figure [j] in the seven-atom cluster problem. 
We use natural units and the same parameters as in (for 
which, e.g., k B T/AE A ^B = 0.033 and 7/2^1^71 = 0.012 - 
low friction limit). The rates kA->B, ks-*A, and fc_B-»c corre- 
spond respectively to the rates for the subprocesses Co — > Cf , 
Cf -> Co, and C'f -> Cf identified in §. The values labeled 
"string" were obtained by the noisy string method using (|^), 
(pd[), and (p"o|). The values labeled "exact" were obtained us- 
ing ( |lC)|) and (|l^), by identifying minima and saddle points 
along the transition path, computing the corresponding en- 
ergy barrier AE, and evaluating all the eigenvalues of the 
Hessian at the minima and the saddle points from the Hes- 
sian itself. 
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FIG. 3: Two of the paths [(a) and (b)] followed by the magnetization vector m during a switching. The pictures show the 
succession minimum - saddle - ... - saddle - minimum. The out-of-plane component of m is very small (less than 10 -2 ) during 
the switching and we only plot its in-plane component with color coding: blue = right, red = left, yellow = up, green = down. 
For both paths, we used N = 200 discretization points along the string. 




FIG. 4: The magnetic energy along the two paths [(a) and 
(b)] shown in figure 



trinsic parametrization of the string leads to a simple 
and efficient algorithm for the numerical solution of its 
evolution equation, and permits to sample the configura- 
tion space in regions that otherwise would be practically 
inaccessible by standard Monte-Carlo methods. 
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